---
title: "No One Mourns the Wicked: The Limits of Partisan Hostility Persisting through Tragedy"
author: "Wayde Z.C. Marsh"
date: "10/31/24"
output:
  html_document:
    toc: true
    code_folding: hide
  latex_engine: xelatex
  number_sections: yes
mainfont: Helvetica
---

# Setting up the Data {.tabset}

## Loading packages and data

```{r setup, message = FALSE, warning = FALSE, results ='hide'}

##### Packages #####

library(descr)
library(ggplot2)
library(dplyr)
library(foreign)
library(plyr)
library(ggpubr)
library(ggeffects)
library(tidyverse)
library(haven)
library(modelsummary)
library(Rmisc)
library(xtable)
options("modelsummary_format_numeric_latex" = "mathmode")


# setwd(dirname(rstudioapi::getSourceEditorContext()$path)) 
# knitr::opts_knit$set(root.dir = '~/')

```

# Raw Data

```{r raw data}

# load Lucid data

df <- read.csv('CES21_pilot1.csv')

# clean Lucid Data

df$Ind.Leaners <- as.numeric(df$Ind.Leaners)
df$Dem.Party.Strength <- as.numeric(df$Dem.Party.Strength)
df$Rep.Party.Strength <- as.numeric(df$Rep.Party.Strength)

df$pid <- ifelse((df$Party.ID.1 == 1 & df$Dem.Party.Strength == 1), 1,
          ifelse((df$Party.ID.1 == 1 & df$Dem.Party.Strength == 2), 2,
          ifelse((df$Party.ID.1 == 3 & df$Ind.Leaners == 2) | 
                 (df$Party.ID.1 == 4 & df$Ind.Leaners == 2) |
                 (df$Party.ID.1 == 5 & df$Ind.Leaners == 2), 3,
          ifelse((df$Party.ID.1 == 3 & df$Ind.Leaners == 3) | 
                 (df$Party.ID.1 == 4 & df$Ind.Leaners == 3) |
                 (df$Party.ID.1 == 5 & df$Ind.Leaners == 3) |
                 (df$Party.ID.1 == 3 & df$Ind.Leaners == 4) | 
                 (df$Party.ID.1 == 4 & df$Ind.Leaners == 4) |
                 (df$Party.ID.1 == 5 & df$Ind.Leaners == 4), 4,
          ifelse((df$Party.ID.1 == 3 & df$Ind.Leaners == 1) | 
                 (df$Party.ID.1 == 4 & df$Ind.Leaners == 1) |
                 (df$Party.ID.1 == 5 & df$Ind.Leaners == 1), 5,
          ifelse((df$Party.ID.1 == 2 & df$Rep.Party.Strength == 2), 6,
          ifelse((df$Party.ID.1 == 2 & df$Rep.Party.Strength == 1), 7, NA)))))))

df$gop <- ifelse(df$pid > 4, 1,
          ifelse(df$pid < 4, 0, NA))

## Creating treatment variable

df$tr2_ct <- ifelse(is.na(df$study2_control_time_Page.Submit), 0, 1)
df$tr2_bl <- ifelse(is.na(df$study2_black_time_Page.Submit), 0, 1)
df$tr2_lx <- ifelse(is.na(df$study2_latino_time_Page.Submit), 0, 1)
df$tr2_wh <- ifelse(is.na(df$study2_white_time_Page.Submit), 0, 1)

df$tr2 <- ifelse(df$tr2_ct == 1, 1,
          ifelse(df$tr2_bl == 1, 2,
          ifelse(df$tr2_lx == 1, 3,
          ifelse(df$tr2_wh == 1, 4, NA))))

## Creating treatment variable

df$tr1_ct <- ifelse(is.na(df$study1_control_time_Page.Submit), 0, 1)
df$tr1_anti <- ifelse(is.na(df$study1_antivax_time_Page.Submit), 0, 1)
df$tr1_dem <- ifelse(is.na(df$study1_dem_time_Page.Submit), 0, 1)
df$tr1_rep <- ifelse(is.na(df$study1_rep_time_Page.Submit), 0, 1)
df$tr1_antidem <- ifelse(is.na(df$study1_antidem_time_Page.Submit), 0, 1)
df$tr1_antirep <- ifelse(is.na(df$study1_antirep_time_Page.Submit), 0, 1)

df$tr1 <- ifelse(df$tr1_ct == 1, 1,
          ifelse(df$tr1_anti == 1, 2,
          ifelse(df$tr1_dem == 1, 3,
          ifelse(df$tr1_rep == 1, 4,
          ifelse(df$tr1_antidem == 1, 5,
          ifelse(df$tr1_antirep == 1, 6, NA))))))

df$CaseID <- 1:nrow(df)

## Combine DVs

df$study2_dv1 <- rowSums(cbind(df$study2_dv_control_1,
                                df$study2_dv_black_1,
                                df$study2_dv_latino_1,
                                df$study2_dv_white_1), na.rm = T)
df$study2_dv2 <- rowSums(cbind(df$study2_dv_control_2,
                                df$study2_dv_black_2,
                                df$study2_dv_latino_2,
                                df$study2_dv_white_2), na.rm = T)
df$study2_dv3 <- rowSums(cbind(df$study2_dv_control_3,
                                df$study2_dv_black_3,
                                df$study2_dv_latino_3,
                                df$study2_dv_white_3), na.rm = T)
df$study2_dv4 <- rowSums(cbind(df$study2_dv_control_4,
                                df$study2_dv_black_4,
                                df$study2_dv_latino_4,
                                df$study2_dv_white_4), na.rm = T)
df$study2_dv_sym <- rowSums(cbind(df$study2_dv2_control_1,
                                df$study2_dv2_black_1,
                                df$study2_dv2_latino_1,
                                df$study2_dv2_white_1), na.rm = T)

# load CES data

df2 <- read_sav('CCES21_NDB_OUTPUT.sav')

# clean CES data

df2$trgrp <- df2$UND452

df2$study1_dv_1 <- df2$UND453
df2$study1_dv_2 <- df2$UND454
df2$study1_dv_3 <- df2$UND455
df2$study1_dv2_1 <- df2$UND456

df2$mandate <- df2$UND448
df2$covidconcern <- df2$UND449
df2$vaxstatus <- df2$CC21_306

df2$pid <- df2$pid7
df2$pid <- ifelse(df2$pid > 7, NA, df2$pid)
df2$gop <- ifelse(df2$pid > 4, 1,
                 ifelse(df2$pid < 4, 0, NA))

df2$coviddeath <- ifelse((df2$CC21_305b_1 == 1 | df2$CC21_305b_2 == 1) |
                           df2$CC21_305b_3 == 1, 1, 2)
df2$coviddeath <- ifelse(is.na(df2$coviddeath), 2, df2$coviddeath)

df$black <- ifelse(df$ethnicity == 2, 1, 0)
df2$black <- ifelse(df2$race == 2, 1, 0)
df$hispanic <- ifelse(df$hispanic == 2, 1, 0)
df2$hispanic <- ifelse(df2$hispanic == 1, 1, 0)
df$female <- ifelse(df$gender == 2, 1, 0)
df2$female <- ifelse(df2$gender4 == 2, 1, 0)

df2$age <- 2021 - df2$birthyr

df$education <- ifelse(df$education < 0, NA,
                ifelse(df$education == 3, 2,
                ifelse((df$education == 4 | df$education == 5), 3,
                ifelse(df$education == 6, 4,
                ifelse(df$education == 7, 5, 
                ifelse(df$education == 8, 6, df$education))))))

df2$covidconcern <- df2$UND449
df2$vaxstatus <- df2$CC21_306

df <- df[, c('CaseID', 'age', 'female', 'black', 'hispanic', 'education',
             'tr1', 'study1_dv_1', 'study1_dv_2', 'study1_dv2_1', 'study1_dv_3',
             'pid', 'gop', 'coviddeath', 'vaxstatus', 'covidconcern_1', 'covidvaccine_1')]

df2 <- df2[, c('caseid', 'age', 'female', 'black', 'hispanic', 'educ',
             'trgrp', 'study1_dv_1', 'study1_dv_2', 'study1_dv2_1', 'study1_dv_3',
             'pid', 'gop', 'teamweight', 'coviddeath', 'vaxstatus', 'covidconcern', 'mandate')]

prop <- function(x, na.rm = FALSE) (x /100) # create function to turn 0-100 scales to 0-1

df <- df %>%
  mutate(across(c('study1_dv_1', 'study1_dv_2', 'study1_dv2_1', 'study1_dv_3',
                  'covidconcern_1', 'covidvaccine_1'), prop)) # convert 0-100 to 0-1 in study 1
df2 <- df2 %>%
  mutate(across(c('study1_dv_1', 'study1_dv_2', 'study1_dv2_1', 'study1_dv_3',
                  'covidconcern', 'mandate'), prop)) # convert 0-100 to 0-1 in study 2

# Create treatment variables

df$tr_dem <- ifelse(df$tr1 == 3 | df$tr1 == 5, 1, 0)
df$tr_rep <- ifelse(df$tr1 == 4 | df$tr1 == 6, 1, 0)
df$tr_pid <- ifelse(df$tr1 == 3 | df$tr1 == 5, 1, 
             ifelse(df$tr1 == 4 | df$tr1 == 6, 2, 0))
df$tr_vax <- ifelse(df$tr1 == 2 | df$tr1 > 4, 1, 0)

df2$tr_dem <- ifelse(df2$trgrp == 3 | df2$trgrp == 5, 1, 0)
df2$tr_rep <- ifelse(df2$trgrp == 4 | df2$trgrp == 6, 1, 0)
df2$tr_pid <- ifelse(df2$trgrp == 3 | df2$trgrp == 5, 1, 
             ifelse(df2$trgrp == 4 | df2$trgrp == 6, 2, 0))
df2$tr_vax <- ifelse(df2$trgrp == 2 | df2$trgrp > 4, 1, 0)

df$coviddeath <- ifelse(df$coviddeath == 2, 0, df$coviddeath) # simplifies covid deat and vaccination status variables
df2$coviddeath <- ifelse(df2$coviddeath == 2, 0, df2$coviddeath)
df2$vaxstatus <- ifelse(df2$vaxstatus == 2, 0, df2$vaxstatus)
df$vaxstatus2 <- ifelse(df$vaxstatus == 3, 0,
                 ifelse(df$vaxstatus > 4, NA, 1))

```

# Analysis

## Modeling Effects

```{r Modeling Effects}
# scale <- function(x){x/100}
# vars <- c('study1_dv_1', 'study1_dv_2', 'study1_dv2_1', 'study1_dv_3')
# df <- df %>% 
#   mutate_at(vars, scale)
# df2 <- df2 %>% 
#   mutate_at(vars, scale)

## Main Analysis (3x2 design)
# models for table 2
# study I
mod1_1 <- lm(study1_dv_1 ~ factor(tr_pid)*factor(tr_vax) + pid, data = df) 
summary(mod1_1)

mod1_dem_1 <- lm(study1_dv_1 ~ factor(tr_pid)*factor(tr_vax), data = subset(df, gop == 0)) 
summary(mod1_dem_1)

mod1_gop_1 <- lm(study1_dv_1 ~ factor(tr_pid)*factor(tr_vax), data = subset(df, gop == 1)) 
summary(mod1_gop_1)

mod1_ind_1 <- lm(study1_dv_1 ~ factor(tr_pid)*factor(tr_vax), data = subset(df, pid == 4)) 
summary(mod1_ind_1)

mod1_ind2_1 <- lm(study1_dv_1 ~ factor(tr_pid)*factor(tr_vax), data = subset(df, pid > 2 & pid < 6)) 
summary(mod1_ind2_1)

# creates table 2
modelsummary(models = list('Full Sample'             = mod1_1,
                           'Democrats'               = mod1_dem_1,
                           'Republicans'             = mod1_gop_1,
                           'Independents'            = mod1_ind_1,
                           'Independents w/ Leaners' = mod1_ind2_1),
             output = 'latex',
             stars = c('*' = .05),
             gof_omit = 'AIC|BIC|F|RMSE|Log.Lik.',
             escape = FALSE)
# models for table 4
mod2_1 <- lm(study1_dv2_1 ~ factor(tr_pid)*factor(tr_vax) + pid, data = df) 
summary(mod2_1)

mod2_dem_1 <- lm(study1_dv2_1 ~ factor(tr_pid)*factor(tr_vax), data = subset(df, gop == 0)) 
summary(mod2_dem_1)

mod2_gop_1 <- lm(study1_dv2_1 ~ factor(tr_pid)*factor(tr_vax), data = subset(df, gop == 1)) 
summary(mod2_gop_1)

mod2_ind_1 <- lm(study1_dv2_1 ~ factor(tr_pid)*factor(tr_vax), data = subset(df, pid == 4)) 
summary(mod2_ind_1)

mod2_ind2_1 <- lm(study1_dv2_1 ~ factor(tr_pid)*factor(tr_vax), data = subset(df, pid > 2 & pid < 6)) 
summary(mod2_ind2_1)
# creates table 4
modelsummary(models = list('Full Sample'             = mod2_1,
                           'Democrats'               = mod2_dem_1,
                           'Republicans'             = mod2_gop_1,
                           'Independents'            = mod2_ind_1,
                           'Independents w/ Leaners' = mod2_ind2_1),
             output = 'latex',
             stars = c('*' = .05),
             gof_omit = 'AIC|BIC|F|RMSE|Log.Lik.',
             escape = FALSE)

# study II
# models for table 3
mod1_2 <- lm(study1_dv_1 ~ factor(tr_pid)*factor(tr_vax) + pid, data = df2, weights = teamweight) 
summary(mod1_2)

mod1_dem_2 <- lm(study1_dv_1 ~ factor(tr_pid)*factor(tr_vax), data = subset(df2, gop == 0), weights = teamweight) 
summary(mod1_dem_2)

mod1_gop_2 <- lm(study1_dv_1 ~ factor(tr_pid)*factor(tr_vax), data = subset(df2, gop == 1), weights = teamweight) 
summary(mod1_gop_2)

mod1_ind_2 <- lm(study1_dv_1 ~ factor(tr_pid)*factor(tr_vax), data = subset(df2, pid == 4), weights = teamweight) 
summary(mod1_ind_2)

mod1_ind2_2 <- lm(study1_dv_1 ~ factor(tr_pid)*factor(tr_vax), data = subset(df2, pid > 2 & pid < 6), weights = teamweight) 
summary(mod1_ind2_2)
# creates table 3
modelsummary(models = list('Full Sample'             = mod1_2,
                           'Democrats'               = mod1_dem_2,
                           'Republicans'             = mod1_gop_2,
                           'Independents'            = mod1_ind_2,
                           'Independents w/ Leaners' = mod1_ind2_2),
             output = 'latex',
             stars = c('*' = .05),
             gof_omit = 'AIC|BIC|F|RMSE|Log.Lik.',
             escape = FALSE)
# models for table 5
mod2_2 <- lm(study1_dv2_1 ~ factor(tr_pid)*factor(tr_vax) + pid, data = df2, weights = teamweight) 
summary(mod2_2)

mod2_dem_2 <- lm(study1_dv2_1 ~ factor(tr_pid)*factor(tr_vax), data = subset(df2, gop == 0), weights = teamweight) 
summary(mod2_dem_2)

mod2_gop_2 <- lm(study1_dv2_1 ~ factor(tr_pid)*factor(tr_vax), data = subset(df2, gop == 1), weights = teamweight) 
summary(mod2_gop_2)

mod2_ind_2 <- lm(study1_dv2_1 ~ factor(tr_pid)*factor(tr_vax), data = subset(df2, pid == 4), weights = teamweight) 
summary(mod2_ind_2)

mod2_ind2_2 <- lm(study1_dv2_1 ~ factor(tr_pid)*factor(tr_vax), data = subset(df2, pid > 2 & pid < 6), weights = teamweight) 
summary(mod2_ind2_2)
# creates table 5
modelsummary(models = list('Full Sample'             = mod2_2,
                           'Democrats'               = mod2_dem_2,
                           'Republicans'             = mod2_gop_2,
                           'Independents'            = mod2_ind_2,
                           'Independents w/ Leaners' = mod2_ind2_2),
             output = 'latex',
             stars = c('*' = .05),
             gof_omit = 'AIC|BIC|F|RMSE|Log.Lik.',
             escape = FALSE)

# Vax Status Models
# models for table 6
# study I
mod_vax <- lm(study1_dv_1 ~ factor(tr_pid)*factor(tr_vax) + pid + vaxstatus2 + coviddeath + covidconcern_1 + covidvaccine_1, data = df)
summary(mod_vax)
mod1_1_novax <- lm(study1_dv_1 ~ factor(tr_pid)*factor(tr_vax) + pid, data = subset(df, vaxstatus2 == 0))
summary(mod1_1_novax)
mod1_1_vax <- lm(study1_dv_1 ~ factor(tr_pid)*factor(tr_vax) + pid, data = subset(df, vaxstatus2 == 1))
summary(mod1_1_vax)

mod1_1_nocov <- lm(study1_dv_1 ~ factor(tr_pid)*factor(tr_vax) + pid, data = subset(df, coviddeath == 0))
summary(mod1_1_nocov)
mod1_1_cov <- lm(study1_dv_1 ~ factor(tr_pid)*factor(tr_vax) + pid, data = subset(df, coviddeath == 1))
summary(mod1_1_cov)

# study II
mod_vax2 <- lm(study1_dv_1 ~ factor(tr_pid)*factor(tr_vax) + pid + vaxstatus + coviddeath + covidconcern + mandate, data = df2, weights = teamweight)
summary(mod_vax2)
mod1_1_novax2 <- lm(study1_dv_1 ~ factor(tr_pid)*factor(tr_vax) + pid, data = subset(df2, vaxstatus == 0), weights = teamweight)
summary(mod1_1_novax2)
mod1_1_vax2 <- lm(study1_dv_1 ~ factor(tr_pid)*factor(tr_vax) + pid, data = subset(df2, vaxstatus == 1), weights = teamweight)
summary(mod1_1_vax2)

mod1_1_nocov2 <- lm(study1_dv_1 ~ factor(tr_pid)*factor(tr_vax) + pid, data = subset(df2, coviddeath == 0), weights = teamweight)
summary(mod1_1_nocov2)
mod1_1_cov2 <- lm(study1_dv_1 ~ factor(tr_pid)*factor(tr_vax) + pid, data = subset(df2, coviddeath == 1), weights = teamweight)
summary(mod1_1_cov2)
# creates table 6
modelsummary(models = list('Study I'  = mod_vax,
                           'Study II' = mod_vax2,
                           'Study I'  = mod1_1_novax,
                           'Study II' = mod1_1_novax2,
                           'Study I'  = mod1_1_vax,
                           'Study II' = mod1_1_vax2,
                           'Study I'  = mod1_1_nocov,
                           'Study II' = mod1_1_nocov2,
                           'Study I'  = mod1_1_cov,
                           'Study II' = mod1_1_cov2),
             estimate = '{estimate}{stars} ({std.error})',
             statistic = '[{p.value}]',
             output = 'latex',
             stars = c('*' = .05),
             gof_omit = 'AIC|BIC|F|RMSE|Log.Lik.',
             escape = FALSE)
# models for table 7
# study I
mod2_vax <- lm(study1_dv2_1 ~ factor(tr_pid)*factor(tr_vax) + pid + vaxstatus2 + coviddeath + covidconcern_1 + covidvaccine_1, data = df)
summary(mod2_vax)
mod2_1_novax <- lm(study1_dv2_1 ~ factor(tr_pid)*factor(tr_vax) + pid, data = subset(df, vaxstatus2 == 0))
summary(mod1_1_novax)
mod2_1_vax <- lm(study1_dv2_1 ~ factor(tr_pid)*factor(tr_vax) + pid, data = subset(df, vaxstatus2 == 1))
summary(mod1_1_vax)

mod2_1_nocov <- lm(study1_dv2_1 ~ factor(tr_pid)*factor(tr_vax) + pid, data = subset(df, coviddeath == 0))
summary(mod1_1_nocov)
mod2_1_cov <- lm(study1_dv2_1 ~ factor(tr_pid)*factor(tr_vax) + pid, data = subset(df, coviddeath == 1))
summary(mod2_1_cov)

# study II
mod2_vax2 <- lm(study1_dv2_1 ~ factor(tr_pid)*factor(tr_vax) + pid + vaxstatus + coviddeath + covidconcern + mandate, data = df2, weights = teamweight)
summary(mod2_vax2)
mod2_1_novax2 <- lm(study1_dv2_1 ~ factor(tr_pid)*factor(tr_vax) + pid, data = subset(df2, vaxstatus == 0), weights = teamweight)
summary(mod2_1_novax2)
mod2_1_vax2 <- lm(study1_dv2_1 ~ factor(tr_pid)*factor(tr_vax) + pid, data = subset(df2, vaxstatus == 1), weights = teamweight)
summary(mod2_1_vax2)

mod2_1_nocov2 <- lm(study1_dv2_1 ~ factor(tr_pid)*factor(tr_vax) + pid, data = subset(df2, coviddeath == 0), weights = teamweight)
summary(mod2_1_nocov2)
mod2_1_cov2 <- lm(study1_dv2_1 ~ factor(tr_pid)*factor(tr_vax) + pid, data = subset(df2, coviddeath == 1), weights = teamweight)
summary(mod2_1_cov2)
# creates table 7
modelsummary(models = list('Study I'  = mod2_vax,
                           'Study II' = mod2_vax2,
                           'Study I'  = mod2_1_novax,
                           'Study II' = mod2_1_novax2,
                           'Study I'  = mod2_1_vax,
                           'Study II' = mod2_1_vax2,
                           'Study I'  = mod2_1_nocov,
                           'Study II' = mod2_1_nocov2,
                           'Study I'  = mod2_1_cov,
                           'Study II' = mod2_1_cov2),
             fmt = fmt_decimal(3, 16),
             estimate = '{estimate}{stars} ({std.error})',
             statistic = '[{p.value}]',
             output = 'latex',
             stars = c('*' = .05),
             gof_omit = 'AIC|BIC|F|RMSE|Log.Lik.',
             escape = FALSE)

### Supplementary Analysis

# blaming victim

m1_1 <- lm(study1_dv_1 ~ factor(tr1), data = df) # 
summary(m1_1)

m1_dem_1 <- lm(study1_dv_1 ~ factor(tr1), data = subset(df, gop == 0)) # 
summary(m1_dem_1)

m1_gop_1 <- lm(study1_dv_1 ~ factor(tr1), data = subset(df, gop == 1)) # 
summary(m1_gop_1)

m2_1 <- lm(study1_dv_1 ~ factor(trgrp), data = df2, weights = teamweight) # 
summary(m2_1)

m2_dem_1 <- lm(study1_dv_1 ~ factor(trgrp), data = subset(df2, gop == 0), weights = teamweight) # 
summary(m2_dem_1)

m2_gop_1 <- lm(study1_dv_1 ~ factor(trgrp), data = subset(df2, gop == 1), weights = teamweight) # 
summary(m2_gop_1)
# creates table F.1 
modelsummary(models = list('Full Sample'  = m1_1,
                           'Democrats' = m1_dem_1,
                           'Republicans' = m1_gop_1),
             output = 'latex',
             stars = c('*' = .05),
             gof_omit = 'AIC|BIC|F|RMSE|Log.Lik.',
             escape = FALSE)
#creates table F.2
modelsummary(models = list('Full Sample'  = m2_1,
                           'Democrats' = m2_dem_1,
                           'Republicans' = m2_gop_1),
             output = 'latex',
             stars = c('*' = .05),
             gof_omit = 'AIC|BIC|F|RMSE|Log.Lik.',
             escape = FALSE)

# blaming govt

m1_2 <- lm(study1_dv_2 ~ factor(tr1), data = df) # 
summary(m1_2)

m1_dem_2 <- lm(study1_dv_2 ~ factor(tr1), data = subset(df, gop == 0)) # 
summary(m1_dem_2)

m1_gop_2 <- lm(study1_dv_2 ~ factor(tr1), data = subset(df, gop == 1)) # 
summary(m1_gop_2)

m2_2 <- lm(study1_dv_2 ~ factor(trgrp), data = df2) # 
summary(m2_2)

m2_dem_2 <- lm(study1_dv_2 ~ factor(trgrp), data = subset(df2, gop == 0)) # 
summary(m2_dem_2)

m2_gop_2 <- lm(study1_dv_2 ~ factor(trgrp), data = subset(df2, gop == 1)) # 
summary(m2_gop_2)

# blaming activists

m1_3 <- lm(study1_dv_3 ~ factor(tr1), data = df) # 
summary(m1_3)

m1_dem_3 <- lm(study1_dv_3 ~ factor(tr1), data = subset(df, gop == 0)) # 
summary(m1_dem_3)

m1_gop_3 <- lm(study1_dv_3 ~ factor(tr1), data = subset(df, gop == 1)) # 
summary(m1_gop_3)

m2_3 <- lm(study1_dv_3 ~ factor(trgrp), data = df2) # 
summary(m2_3)

m2_dem_3 <- lm(study1_dv_3 ~ factor(trgrp), data = subset(df2, gop == 0)) # 
summary(m2_dem_3)

m2_gop_3 <- lm(study1_dv_3 ~ factor(trgrp), data = subset(df2, gop == 1)) # 
summary(m2_gop_3)

# sympathy

m1_4 <- lm(study1_dv2_1 ~ factor(tr1), data = df) # 
summary(m1_4)

m1_dem_4 <- lm(study1_dv2_1 ~ factor(tr1), data = subset(df, gop == 0)) # 
summary(m1_dem_4)

m1_gop_4 <- lm(study1_dv2_1 ~ factor(tr1), data = subset(df, gop == 1)) # 
summary(m1_gop_4)

m2_4 <- lm(study1_dv2_1 ~ factor(trgrp), data = df2, weights = teamweight) # 
summary(m2_4)

m2_dem_4 <- lm(study1_dv2_1 ~ factor(trgrp), data = subset(df2, gop == 0), weights = teamweight) # 
summary(m2_dem_4)

m2_gop_4 <- lm(study1_dv2_1 ~ factor(trgrp), data = subset(df2, gop == 1), weights = teamweight) # 
summary(m2_gop_4)
# creates table F.3
modelsummary(models = list('Full Sample'  = m1_4,
                           'Democrats' = m1_dem_4,
                           'Republicans' = m1_gop_4),
             output = 'latex',
             stars = c('*' = .05),
             gof_omit = 'AIC|BIC|F|RMSE|Log.Lik.',
             escape = FALSE)
# creates table F.4
modelsummary(models = list('Full Sample'  = m2_4,
                           'Democrats' = m2_dem_4,
                           'Republicans' = m2_gop_4),
             output = 'latex',
             stars = c('*' = .05),
             gof_omit = 'AIC|BIC|F|RMSE|Log.Lik.',
             escape = FALSE)


# Testing Group Empathy Theory
# analysis in SI, section G
## Blame

em_death1.1 <- lm(study1_dv_1 ~ factor(tr1), data = subset(df, coviddeath == 1))
summary(em_death1.1)

em_nodeath1.1 <- lm(study1_dv_1 ~ factor(tr1), data = subset(df, coviddeath == 0))
summary(em_nodeath1.1)

em_death2.1 <- lm(study1_dv_1 ~ factor(trgrp), data = subset(df2, coviddeath == 1), weights = teamweight)
summary(em_death2.1)

em_nodeath2.1 <- lm(study1_dv_1 ~ factor(trgrp), data = subset(df2, coviddeath == 0), weights = teamweight)
summary(em_nodeath2.1)
# table G.1
modelsummary(models = list('COVID Death'  = em_death1.1,
                           'No COVID Death' = em_nodeath1.1,
                           'COVID Death'  = em_death2.1,
                           'No COVID Death' = em_nodeath2.1),
             output = 'latex',
             stars = c('*' = .05),
             gof_omit = 'AIC|BIC|F|RMSE|Log.Lik.',
             escape = FALSE)

## Sympathy

em_death1.2 <- lm(study1_dv2_1 ~ factor(tr1), data = subset(df, coviddeath == 1))
summary(em_death1.2)

em_nodeath1.2 <- lm(study1_dv2_1 ~ factor(tr1), data = subset(df, coviddeath == 0))
summary(em_nodeath1.2)

em_death2.2 <- lm(study1_dv2_1 ~ factor(trgrp), data = subset(df2, coviddeath == 1))
summary(em_death2.2)

em_nodeath2.2 <- lm(study1_dv2_1 ~ factor(trgrp), data = subset(df2, coviddeath == 0))
summary(em_nodeath2.2)
# table G.2
modelsummary(models = list('COVID Death'  = em_death1.2,
                           'No COVID Death' = em_nodeath1.2,
                           'COVID Death'  = em_death2.2,
                           'No COVID Death' = em_nodeath2.2),
             output = 'latex',
             stars = c('*' = .05),
             gof_omit = 'AIC|BIC|F|RMSE|Log.Lik.',
             escape = FALSE)

## Testing by Vax Status

# df$vaxstatus2 <- ifelse((df$vaxstatus == 1 | df$vaxstatus == 2), 1,
#                  ifelse(df$vaxstatus == 3, 0, 2))

vac1.1_vax <- lm(study1_dv_1 ~ factor(tr1), data = subset(df, vaxstatus2 == 1))
summary(vac1.1_vax)

vac1.1_novax <- lm(study1_dv_1 ~ factor(tr1), data = subset(df, vaxstatus2 == 0))
summary(vac1.1_novax)

vac1.2_vax <- lm(study1_dv2_1 ~ factor(tr1), data = subset(df, vaxstatus2 == 1))
summary(vac1.2_vax)

vac1.2_novax <- lm(study1_dv2_1 ~ factor(tr1), data = subset(df, vaxstatus2 == 0))
summary(vac1.2_novax)
# table G.3
modelsummary(models = list('Vaccinated'  = vac1.1_vax,
                           'Not Vaccinated' = vac1.1_novax,
                           'Vaccinated'  = vac1.2_vax,
                           'Not Vaccinated' = vac1.2_novax),
             output = 'latex',
             stars = c('*' = .05),
             gof_omit = 'AIC|BIC|F|RMSE|Log.Lik.',
             escape = FALSE)

```

## DIMs (Tukey Multiple pairwise-comparisons)

```{r}

aov_b1r <- aov(study1_dv_1 ~ as.factor(tr1), data = subset(df, gop == 1))
dim_b1r <- TukeyHSD(aov_b1r)
aov_b2r <- aov(study1_dv_1 ~ as.factor(trgrp), data = subset(df2, gop == 1))
dim_b2r <- TukeyHSD(aov_b2r)

aov_s1r <- aov(study1_dv_2 ~ as.factor(tr1), data = subset(df, gop == 1))
dim_s1r <- TukeyHSD(aov_s1r)
aov_s2r <- aov(study1_dv_2 ~ as.factor(trgrp), data = subset(df2, gop == 1))
dim_s2r <- TukeyHSD(aov_s2r)


aov_b1d <- aov(study1_dv_1 ~ as.factor(tr1), data = subset(df, gop == 0))
dim_b1d <- TukeyHSD(aov_b1d)
aov_b2d <- aov(study1_dv_1 ~ as.factor(trgrp), data = subset(df2, gop == 0))
dim_b2d <- TukeyHSD(aov_b2d)

aov_s1d <- aov(study1_dv_2 ~ as.factor(tr1), data = subset(df, gop == 0))
dim_s1d <- TukeyHSD(aov_s1d)
aov_s2d <- aov(study1_dv_2 ~ as.factor(trgrp), data = subset(df2, gop == 0))
dim_s2d <- TukeyHSD(aov_s2d)

dim_b1r <- tidy(dim_b1r)
dim_b2r <- tidy(dim_b2r)
dim_b1d <- tidy(dim_b1d)
dim_b2d <- tidy(dim_b2d)

dim_s1r <- tidy(dim_s1r)
dim_s2r <- tidy(dim_s2r)
dim_s1d <- tidy(dim_s1d)
dim_s2d <- tidy(dim_s2d)


library(broom) # creates tables E.1 through E.8
print.xtable(xtable(dim_b1r[c(2, 4, 7)]), include.rownames = F)
print.xtable(xtable(dim_b1d[c(2, 4, 7)]), include.rownames = F)
print.xtable(xtable(dim_b2r[c(2, 4, 7)]), include.rownames = F)
print.xtable(xtable(dim_b2d[c(2, 4, 7)]), include.rownames = F)

print.xtable(xtable(dim_s1r[c(2, 4, 7)]), include.rownames = F)
print.xtable(xtable(dim_s1d[c(2, 4, 7)]), include.rownames = F)
print.xtable(xtable(dim_s2r[c(2, 4, 7)]), include.rownames = F)
print.xtable(xtable(dim_s2d[c(2, 4, 7)]), include.rownames = F)

```

# Summary Statistics and Diagnostics

## Survey Stats

```{r}

round(colMeans(df, na.rm = T), 2) # to create table B.3--vaxstatus2 is the proper vaccine status variable to use
format(round(colMeans(df2, na.rm = T), 3), scientific = F) # to create table B.3

# Descriptives of DVs

hist.df <- as.data.frame(rbind(df[, c('study1_dv_1', 'study1_dv2_1')], df2[, c('study1_dv_1', 'study1_dv2_1')]))
hist.df$study <- c(rep('Study I', 1444), rep('Study II', 1000))
names(hist.df) <- c('Blame', 'Sympathy', 'Study')


dv_hist1 <- ggplot(hist.df, aes(x = Blame)) + # creates figure B.6
  geom_histogram(aes(y = ..density..), binwidth = 0.01) +
  scale_y_continuous(labels = scales::percent_format()) +
  facet_wrap(~Study) +
  theme(plot.title = element_text(hjust = 0.5)) +
  theme(legend.position="bottom") +
  theme(text=element_text(family="serif"))

dv_hist2 <- ggplot(hist.df, aes(x = Sympathy)) + # creates figure B.7
  geom_histogram(aes(y = ..density..), binwidth = 0.01) +
  scale_y_continuous(labels = scales::percent_format()) +
  facet_wrap(~Study) +
  theme(plot.title = element_text(hjust = 0.5)) +
  theme(legend.position="bottom") +
  theme(text=element_text(family="serif"))

# ggsave('dv_hist1.png', dv_hist1, 'png', '~/Dropbox/Wayde/Research/Apps/Overleaf/CovidPolar_JEPS',
#     width = 190, height = 127, units = 'mm')
# ggsave('dv_hist2.png', dv_hist2, 'png', '~/Dropbox/Wayde/Research/Apps/Overleaf/CovidPolar_JEPS',
#     width = 190, height = 127, units = 'mm')

```

## Means of DVs

```{r}

m1_blame <- aggregate(study1_dv_1 ~ tr1, data = df, FUN = mean)
m1_symp <- aggregate(study1_dv_2 ~ tr1, data = df, FUN = mean)
s1_blame <- aggregate(study1_dv_1 ~ tr1, data = df, FUN = sd)
s1_symp <- aggregate(study1_dv_2 ~ tr1, data = df, FUN = sd)

m2_blame <- aggregate(study1_dv_1 ~ trgrp, data = df2, FUN = mean)
m2_symp <- aggregate(study1_dv_2 ~ trgrp, data = df2, FUN = mean)
s2_blame <- aggregate(study1_dv_1 ~ trgrp, data = df2, FUN = sd)
s2_symp <- aggregate(study1_dv_2 ~ trgrp, data = df2, FUN = sd)

names(m1_blame) <- c('trgrp', 'mean')
names(m1_symp) <- c('trgrp', 'mean')
names(m2_blame) <- c('trgrp', 'mean')
names(m2_symp) <- c('trgrp', 'mean')
names(s1_blame) <- c('trgrp', 'sd')
names(s1_symp) <- c('trgrp', 'sd')
names(s2_blame) <- c('trgrp', 'sd')
names(s2_symp) <- c('trgrp', 'sd')

names_vec <- rbind(m1_blame[1], m2_blame[1],
                m1_symp[1], m2_symp[1])
means_vec <- rbind(m1_blame[2], m2_blame[2],
                m1_symp[2], m2_symp[2])
sd_vec <- rbind(s1_blame[2], s2_blame[2],
                s1_symp[2], s2_symp[2])
sum_df <- cbind(names_vec, means_vec, sd_vec)

sum_df <- sum_df %>%
  mutate(trgrp = recode(trgrp, '1' = 'Control', '2' = 'Anti-Vax', '3' = 'Democrat', '4' = 'Republican',
                        '5' = 'Anti-Vax Democrat', '6' = 'Anti-Vax Republican'))
sum_df$dv <- c(rep('Blame DV', 12), rep('Sympathy DV', 12))
sum_df$study <- c(rep('Study I', 6), rep('Study II', 6),
                  rep('Study I', 6), rep('Study II', 6))

library(xtable) # tables B.1 and B.2
print.xtable(xtable(sum_df[sum_df$dv == 'Blame DV', c(1:3, 5)], type = 'latex'), include.rownames = F)
print.xtable(xtable(sum_df[sum_df$dv == 'Sympathy DV', c(1:3, 5)], type = 'latex'), include.rownames = F)

```

## Balance Analysis

```{r}

df_complete <- na.omit(df)
df2_complete <- na.omit(df2)

model1_age <- aov(age ~ as.factor(tr1), data = df_complete)
model1_edu <- aov(education ~ as.factor(tr1), data = df_complete)
model1_sex <- aov(female ~ as.factor(tr1), data = df_complete)
model1_blk <- aov(black ~ as.factor(tr1), data = df_complete)
model1_hsp <- aov(hispanic ~ as.factor(tr1), data = df_complete)
model1_pid <- aov(pid ~ as.factor(tr1), data = df_complete)

model2_age <- aov(age ~ as.factor(trgrp), data = df2_complete)
model2_edu <- aov(educ ~ as.factor(trgrp), data = df2_complete)
model2_sex <- aov(female ~ as.factor(trgrp), data = df2_complete)
model2_blk <- aov(black ~ as.factor(trgrp), data = df2_complete)
model2_hsp <- aov(hispanic ~ as.factor(trgrp), data = df2_complete)
model2_pid <- aov(pid ~ as.factor(trgrp), data = df2_complete)

gg_tukey<-function(Tukey){
  A<-require("tidyverse")
  if(A==TRUE){
    library(tidyverse)
  } else {
    install.packages("tidyverse")
    library(tidyverse)
  }
  B<-as.data.frame(Tukey[1])
  colnames(B)[2:3]<-c("min",
                      "max")
  C<-data.frame(id=row.names(B),
                min=B$min,
                max=B$max)
  D<-C%>%
    ggplot(aes(id))+
    geom_errorbar(aes(ymin=min,
                     ymax=max),
                 width = 0.2)+
    geom_hline(yintercept=0,
               color="red")+
    labs(x=NULL)+
    coord_flip()+
    theme(text=element_text(family="serif"),
            title=element_text(color="black",size=15),
            axis.text = element_text(color="black",size=10),
            axis.title = element_text(color="black",size=10),
            panel.grid=element_line(color="grey75"),
            axis.line=element_blank(),
            plot.background=element_rect(fill="white",color="white"),
            panel.background=element_rect(fill="white"),
            panel.border = element_rect(colour = "black", fill = NA,size=0.59),
            legend.key= element_rect(color="white",fill="white")
          )
  return(D)
}

tk1_1 <- gg_tukey(TukeyHSD(model1_age))
tk1_2 <- gg_tukey(TukeyHSD(model1_edu))
tk1_3 <- gg_tukey(TukeyHSD(model1_sex))
tk1_4 <- gg_tukey(TukeyHSD(model1_blk))
tk1_5 <- gg_tukey(TukeyHSD(model1_hsp))
tk1_6 <- gg_tukey(TukeyHSD(model1_pid))

tk2_1 <- gg_tukey(TukeyHSD(model2_age))
tk2_2 <- gg_tukey(TukeyHSD(model2_edu))
tk2_3 <- gg_tukey(TukeyHSD(model2_sex))
tk2_4 <- gg_tukey(TukeyHSD(model2_blk))
tk2_5 <- gg_tukey(TukeyHSD(model2_hsp))
tk2_6 <- gg_tukey(TukeyHSD(model2_pid))

# tk1_1 <- plot(TukeyHSD(model1_age, conf.level = 0.95), las = 2)
# tk1_2 <- plot(TukeyHSD(model1_edu, conf.level = 0.95), las = 2)
# tk1_3 <- plot(TukeyHSD(model1_sex, conf.level = 0.95), las = 2)
# tk1_4 <- plot(TukeyHSD(model1_blk, conf.level = 0.95), las = 2)
# tk1_5 <- plot(TukeyHSD(model1_hsp, conf.level = 0.95), las = 2)
# tk1_6 <- plot(TukeyHSD(model1_pid, conf.level = 0.95), las = 2)
# 
# tk2_1 <- plot(TukeyHSD(model2_age, conf.level = 0.95), las = 2)
# tk2_2 <- plot(TukeyHSD(model2_edu, conf.level = 0.95), las = 2)
# tk2_3 <- plot(TukeyHSD(model2_sex, conf.level = 0.95), las = 2)
# tk2_4 <- plot(TukeyHSD(model2_blk, conf.level = 0.95), las = 2)
# tk2_5 <- plot(TukeyHSD(model2_hsp, conf.level = 0.95), las = 2)
# tk2_6 <- plot(TukeyHSD(model2_pid, conf.level = 0.95), las = 2)

tk1 <- ggarrange(tk1_1, tk1_2, tk1_3, tk1_4, tk1_5, tk1_6,# creates figure B.8
          labels = c('Age', 'Education', 'Sex', 'Black', 'Hispanic', 'Party ID'),
          ncol = 2, nrow = 3)

tk2 <- ggarrange(tk2_1, tk2_2, tk2_3, tk2_4, tk2_5, tk2_6,# creates figure B.9
          labels = c('Age', 'Education', 'Sex', 'Black', 'Hispanic', 'Party ID'),
          ncol = 2, nrow = 3)


```

# Additional/Secondary Figures and Analysis

## Figures

```{r Figures}

# blaming victim

coef1 <- c(m1_dem_1$coefficients[2:6], m1_gop_1$coefficients[2:6],
           m2_dem_1$coefficients[2:6], m2_gop_1$coefficients[2:6])
se1 <- sqrt(diag(vcov(m1_dem_1))); se2 <- sqrt(diag(vcov(m1_gop_1))); se3 <- sqrt(diag(vcov(m2_dem_1))); se4 <- sqrt(diag(vcov(m2_gop_1)))

hi1 <- c(coef1[1]+(1.96*se1[2]), coef1[2]+(1.96*se1[3]),
            coef1[3]+(1.96*se1[4]), coef1[4]+(1.96*se1[5]),
            coef1[5]+(1.96*se1[6]),
         coef1[6]+(1.96*se2[2]), coef1[7]+(1.96*se2[3]),
            coef1[8]+(1.96*se2[4]), coef1[9]+(1.96*se2[5]),
            coef1[10]+(1.96*se2[6]),
         coef1[11]+(1.96*se3[2]), coef1[12]+(1.96*se3[3]),
            coef1[13]+(1.96*se3[4]), coef1[14]+(1.96*se3[5]),
            coef1[15]+(1.96*se3[6]),
         coef1[16]+(1.96*se4[2]), coef1[17]+(1.96*se4[3]),
            coef1[18]+(1.96*se4[4]), coef1[19]+(1.96*se4[5]),
            coef1[20]+(1.96*se4[6]))

lo1 <- c(coef1[1]-(1.96*se1[2]), coef1[2]-(1.96*se1[3]),
            coef1[3]-(1.96*se1[4]), coef1[4]-(1.96*se1[5]),
            coef1[5]-(1.96*se1[6]),
         coef1[6]-(1.96*se2[2]), coef1[7]-(1.96*se2[3]),
            coef1[8]-(1.96*se2[4]), coef1[9]-(1.96*se2[5]),
            coef1[10]-(1.96*se2[6]),
         coef1[11]-(1.96*se3[2]), coef1[12]-(1.96*se3[3]),
            coef1[13]-(1.96*se3[4]), coef1[14]-(1.96*se3[5]),
            coef1[15]-(1.96*se3[6]),
         coef1[16]-(1.96*se4[2]), coef1[17]-(1.96*se4[3]),
            coef1[18]-(1.96*se4[4]), coef1[19]-(1.96*se4[5]),
            coef1[20]-(1.96*se4[6]))

study <- c(rep('Study I', 10), rep('Study II', 10))
party <- rep(c('Democratic', 'Democratic', 'Democratic', 'Democratic', 'Democratic',
               'Republican', 'Republican', 'Republican', 'Republican',
               'Republican'), 2)
tr <- rep(c('eAnti-Vax', 'dDemocrat', 'cRepublican',
            'bAnti-Vax Democrat', 'aAnti-Vax Republican'), 4)

df_1 <- as.data.frame(cbind(coef1, hi1, lo1, party, tr, study))
df_1[, 1:3] <- lapply(df_1[, 1:3], as.numeric)
df_1[, 4:6] <- lapply(df_1[, 4:6], as.factor)
df_1$stupid <- str_c(df_1$party, ', ', df_1$study)


plot1 <- ggplot(df_1, aes(shape = party)) +
  theme_bw() +
  geom_vline(xintercept = 0, linetype = 'dashed') +
  geom_pointrange(aes(x = coef1, y = tr, xmin = lo1,
                         xmax = hi1), size = 1,
                     lwd = 1, position = position_dodge(width = 0.5), fill = 'WHITE') +
  scale_shape_manual(name = 'Partisanship',
                    labels = c('Democrats', 'Republicans'),
                    values = c(15, 1)) +
  scale_y_discrete(name = 'Treatment Condition',
                   labels = c('Anti-Vax Republican', 'Anti-Vax Democrat',
                              'Republican', 'Democrat', 'Anti-Vax')) +
  scale_x_continuous("Estimated Effect on Blame Attributed to Victim") +
  facet_wrap(~study) +
  labs(title = 'Partisans Blame Deceased Anti-Vaxxers and Outpartisans for Own Death',
       caption = '\nNote: Estimated effects were calculated from two separate models (based on partisan identity). \nError bars represent 95/% confidence intervals. Full results table in Appendix E.') +
  theme(axis.text=element_text(size = 10)) + 
  theme(plot.title = element_text(hjust = 0.5)) + 
  theme(plot.caption = element_text(hjust = 0)) +
  theme(text=element_text(size = 10, family="serif"))

plot1_sl <- ggplot(df_1, aes(shape = party)) +
  theme_bw() +
  geom_vline(xintercept = 0, linetype = 'dashed') +
  geom_pointrange(aes(x = coef1, y = tr, xmin = lo1,
                         xmax = hi1), size = 1,
                     lwd = 1, position = position_dodge(width = 0.5), fill = 'WHITE') +
  scale_shape_manual(name = 'Partisanship',
                    labels = c('Democrats', 'Republicans'),
                    values = c(15, 1)) +
  scale_y_discrete(name = 'Treatment Condition',
                   labels = c('Anti-Vax Republican', 'Anti-Vax Democrat',
                              'Republican', 'Democrat', 'Anti-Vax')) +
  scale_x_continuous("Estimated Effect on Blame Attributed to Victim") +
  facet_wrap(~study) +
  labs(title = 'Partisans Blame Deceased Anti-Vaxxers and Outpartisans for Own Death') +
  theme(axis.text=element_text(size = 10)) + 
  theme(plot.title = element_text(hjust = 0.5)) + 
  theme(plot.caption = element_text(hjust = 0)) +
  theme(text=element_text(size = 10, family="serif"))

# sympathy

coef2 <- c(m1_dem_4$coefficients[2:6], m1_gop_4$coefficients[2:6],
           m2_dem_4$coefficients[2:6], m2_gop_4$coefficients[2:6])
se1 <- sqrt(diag(vcov(m1_dem_4))); se2 <- sqrt(diag(vcov(m1_gop_4))); se3 <- sqrt(diag(vcov(m2_dem_4))); se4 <- sqrt(diag(vcov(m2_gop_4)))

hi2 <- c(coef2[1]+(1.96*se1[2]), coef2[2]+(1.96*se1[3]),
            coef2[3]+(1.96*se1[4]), coef2[4]+(1.96*se1[5]),
            coef2[5]+(1.96*se1[6]),
         coef2[6]+(1.96*se2[2]), coef2[7]+(1.96*se2[3]),
            coef2[8]+(1.96*se2[4]), coef2[9]+(1.96*se2[5]),
            coef2[10]+(1.96*se2[6]),
         coef2[11]+(1.96*se3[2]), coef2[12]+(1.96*se3[3]),
            coef2[13]+(1.96*se3[4]), coef2[14]+(1.96*se3[5]),
            coef2[15]+(1.96*se3[6]),
         coef2[16]+(1.96*se4[2]), coef2[17]+(1.96*se4[3]),
            coef2[18]+(1.96*se4[4]), coef2[19]+(1.96*se4[5]),
            coef2[20]+(1.96*se4[6]))

lo2 <- c(coef2[1]-(1.96*se1[2]), coef2[2]-(1.96*se1[3]),
            coef2[3]-(1.96*se1[4]), coef2[4]-(1.96*se1[5]),
            coef2[5]-(1.96*se1[6]),
         coef2[6]-(1.96*se2[2]), coef2[7]-(1.96*se2[3]),
            coef2[8]-(1.96*se2[4]), coef2[9]-(1.96*se2[5]),
            coef2[10]-(1.96*se2[6]),
         coef2[11]-(1.96*se3[2]), coef2[12]-(1.96*se3[3]),
            coef2[13]-(1.96*se3[4]), coef2[14]-(1.96*se3[5]),
            coef2[15]-(1.96*se3[6]),
         coef2[16]-(1.96*se4[2]), coef2[17]-(1.96*se4[3]),
            coef2[18]-(1.96*se4[4]), coef2[19]-(1.96*se4[5]),
            coef2[20]-(1.96*se4[6]))

df_2 <- as.data.frame(cbind(coef2, hi2, lo2, party, tr, study))
df_2[, 1:3] <- lapply(df_2[, 1:3], as.numeric)
df_2[, 4:6] <- lapply(df_2[, 4:6], as.factor)
df_2$stupid <- str_c(df_2$party, ', ', df_2$study)


plot2 <- ggplot(df_2, aes(shape = party)) +
  theme_bw() +
  geom_vline(xintercept = 0, linetype = 'dashed') +
  geom_pointrange(aes(x = coef2, y = tr, xmin = lo2,
                         xmax = hi2), size = 1,
                     lwd = 1, position = position_dodge(width = 0.5), fill = 'WHITE') +
  scale_shape_manual(name = 'Partisanship',
                    labels = c('Democrats', 'Republicans'),
                    values = c(15, 1)) +
  scale_y_discrete(name = 'Treatment Condition',
                   labels = c('Anti-Vax Republican', 'Anti-Vax Democrat',
                              'Republican', 'Democrat', 'Anti-Vax')) +
  scale_x_continuous("Estimated Effect on Sympathy") +
  facet_wrap(~study) +
  labs(title = 'Partisanship Structures Sympathy for Covid Victims',
       caption = '\nNote: Estimated effects were calculated from two separate models (based on partisan identity). \nError bars represent 95/% confidence intervals. Full results table in Appendix E.') +
  theme(axis.text=element_text(size = 10)) + 
  theme(plot.title = element_text(hjust = 0.5)) + 
  theme(plot.caption = element_text(hjust = 0)) +
  theme(text=element_text(size = 10, family="serif"))

plot2_sl <- ggplot(df_2, aes(shape = party)) +
  theme_bw() +
  geom_vline(xintercept = 0, linetype = 'dashed') +
  geom_pointrange(aes(x = coef2, y = tr, xmin = lo2,
                         xmax = hi2), size = 1,
                     lwd = 1, position = position_dodge(width = 0.5), fill = 'WHITE') +
  scale_shape_manual(name = 'Partisanship',
                    labels = c('Democrats', 'Republicans'),
                    values = c(15, 1)) +
  scale_y_discrete(name = 'Treatment Condition',
                   labels = c('Anti-Vax Republican', 'Anti-Vax Democrat',
                              'Republican', 'Democrat', 'Anti-Vax')) +
  scale_x_continuous("Estimated Effect on Sympathy") +
  facet_wrap(~study) +
  labs(title = 'Partisanship Structures Sympathy for Covid Victims') +
  theme(axis.text=element_text(size = 10)) + 
  theme(plot.title = element_text(hjust = 0.5)) + 
  theme(plot.caption = element_text(hjust = 0)) +
  theme(text=element_text(size = 10, family="serif"))

# Testing Group Empathy

## Blame

coef1 <- c(em_death1.1$coefficients[2:6], em_nodeath1.1$coefficients[2:6],
           em_death2.1$coefficients[2:6], em_nodeath2.1$coefficients[2:6])
se1 <- sqrt(diag(vcov(em_death1.1))); se2 <- sqrt(diag(vcov(em_nodeath1.1))); se3 <- sqrt(diag(vcov(em_death2.1))); se4 <- sqrt(diag(vcov(em_nodeath2.1)))

hi1 <- c(coef1[1]+(1.96*se1[2]), coef1[2]+(1.96*se1[3]),
            coef1[3]+(1.96*se1[4]), coef1[4]+(1.96*se1[5]),
            coef1[5]+(1.96*se1[6]),
         coef1[6]+(1.96*se2[2]), coef1[7]+(1.96*se2[3]),
            coef1[8]+(1.96*se2[4]), coef1[9]+(1.96*se2[5]),
            coef1[10]+(1.96*se2[6]),
         coef1[11]+(1.96*se3[2]), coef1[12]+(1.96*se3[3]),
            coef1[13]+(1.96*se3[4]), coef1[14]+(1.96*se3[5]),
            coef1[15]+(1.96*se3[6]),
         coef1[16]+(1.96*se4[2]), coef1[17]+(1.96*se4[3]),
            coef1[18]+(1.96*se4[4]), coef1[19]+(1.96*se4[5]),
            coef1[20]+(1.96*se4[6]))

lo1 <- c(coef1[1]-(1.96*se1[2]), coef1[2]-(1.96*se1[3]),
            coef1[3]-(1.96*se1[4]), coef1[4]-(1.96*se1[5]),
            coef1[5]-(1.96*se1[6]),
         coef1[6]-(1.96*se2[2]), coef1[7]-(1.96*se2[3]),
            coef1[8]-(1.96*se2[4]), coef1[9]-(1.96*se2[5]),
            coef1[10]-(1.96*se2[6]),
         coef1[11]-(1.96*se3[2]), coef1[12]-(1.96*se3[3]),
            coef1[13]-(1.96*se3[4]), coef1[14]-(1.96*se3[5]),
            coef1[15]-(1.96*se3[6]),
         coef1[16]-(1.96*se4[2]), coef1[17]-(1.96*se4[3]),
            coef1[18]-(1.96*se4[4]), coef1[19]-(1.96*se4[5]),
            coef1[20]-(1.96*se4[6]))

study <- c(rep('Study I', 10), rep('Study II', 10))
party <- rep(c('No', 'No', 'No', 'No', 'No',
               'Yes', 'Yes', 'Yes', 'Yes', 'Yes'), 2)
tr <- rep(c('eAnti-Vax', 'dDemocrat', 'cRepublican',
            'bAnti-Vax Democrat', 'aAnti-Vax Republican'), 4)

df_3 <- as.data.frame(cbind(coef1, hi1, lo1, party, tr, study))
df_3[, 1:3] <- lapply(df_3[, 1:3], as.numeric)
df_3[, 4:6] <- lapply(df_3[, 4:6], as.factor)
df_3$stupid <- str_c(df_3$party, ', ', df_3$study)

# creating figure G.1
plot3 <- ggplot(df_3, aes(shape = party)) +
  theme_bw() +
  geom_vline(xintercept = 0, linetype = 'dashed') +
  geom_pointrange(aes(x = coef1, y = tr, xmin = lo1,
                         xmax = hi1), size = 1,
                     lwd = 1, position = position_dodge(width = 0.5), fill = 'WHITE') +
  scale_shape_manual(name = 'Close Covid Death',
                    labels = c('Yes', 'No'),
                    values = c(15, 1)) +
  scale_y_discrete(name = 'Treatment Condition',
                   labels = c('Anti-Vax Republican', 'Anti-Vax Democrat',
                              'Republican', 'Democrat', 'Anti-Vax')) +
  scale_x_continuous("Estimated Effect on Blame Attributed to Victim") +
  facet_wrap(~study) +
  labs(title = 'Group Empathy and Blame Attribution',
       caption = '\nNote: Estimated effects were calculated from two separate models (based on partisan identity). \nError bars represent 95/% confidence intervals. Full results table in Appendix E.') +
  theme(axis.text=element_text(size = 10)) + 
  theme(plot.title = element_text(hjust = 0.5)) + 
  theme(plot.caption = element_text(hjust = 0)) +
  theme(text=element_text(size = 10, family="serif"))

plot3_sl <- ggplot(df_3, aes(shape = party)) +
  theme_bw() +
  geom_vline(xintercept = 0, linetype = 'dashed') +
  geom_pointrange(aes(x = coef1, y = tr, xmin = lo1,
                         xmax = hi1), size = 1,
                     lwd = 1, position = position_dodge(width = 0.5), fill = 'WHITE') +
  scale_shape_manual(name = 'Close Covid Death',
                    labels = c('Yes', 'No'),
                    values = c(15, 1)) +
  scale_y_discrete(name = 'Treatment Condition',
                   labels = c('Anti-Vax Republican', 'Anti-Vax Democrat',
                              'Republican', 'Democrat', 'Anti-Vax')) +
  scale_x_continuous("Estimated Effect on Blame Attributed to Victim") +
  facet_wrap(~study) +
  labs(title = 'Group Empathy and Blame Attribution') +
  theme(axis.text=element_text(size = 10)) + 
  theme(plot.title = element_text(hjust = 0.5)) + 
  theme(plot.caption = element_text(hjust = 0)) +
  theme(text=element_text(size = 10, family="serif"))

## Sympathy

coef1 <- c(em_death1.2$coefficients[2:6], em_nodeath1.2$coefficients[2:6],
           em_death2.2$coefficients[2:6], em_nodeath2.2$coefficients[2:6])
se1 <- sqrt(diag(vcov(em_death1.2))); se2 <- sqrt(diag(vcov(em_nodeath1.2))); se3 <- sqrt(diag(vcov(em_death2.2))); se4 <- sqrt(diag(vcov(em_nodeath2.2)))

hi1 <- c(coef1[1]+(1.96*se1[2]), coef1[2]+(1.96*se1[3]),
            coef1[3]+(1.96*se1[4]), coef1[4]+(1.96*se1[5]),
            coef1[5]+(1.96*se1[6]),
         coef1[6]+(1.96*se2[2]), coef1[7]+(1.96*se2[3]),
            coef1[8]+(1.96*se2[4]), coef1[9]+(1.96*se2[5]),
            coef1[10]+(1.96*se2[6]),
         coef1[11]+(1.96*se3[2]), coef1[12]+(1.96*se3[3]),
            coef1[13]+(1.96*se3[4]), coef1[14]+(1.96*se3[5]),
            coef1[15]+(1.96*se3[6]),
         coef1[16]+(1.96*se4[2]), coef1[17]+(1.96*se4[3]),
            coef1[18]+(1.96*se4[4]), coef1[19]+(1.96*se4[5]),
            coef1[20]+(1.96*se4[6]))

lo1 <- c(coef1[1]-(1.96*se1[2]), coef1[2]-(1.96*se1[3]),
            coef1[3]-(1.96*se1[4]), coef1[4]-(1.96*se1[5]),
            coef1[5]-(1.96*se1[6]),
         coef1[6]-(1.96*se2[2]), coef1[7]-(1.96*se2[3]),
            coef1[8]-(1.96*se2[4]), coef1[9]-(1.96*se2[5]),
            coef1[10]-(1.96*se2[6]),
         coef1[11]-(1.96*se3[2]), coef1[12]-(1.96*se3[3]),
            coef1[13]-(1.96*se3[4]), coef1[14]-(1.96*se3[5]),
            coef1[15]-(1.96*se3[6]),
         coef1[16]-(1.96*se4[2]), coef1[17]-(1.96*se4[3]),
            coef1[18]-(1.96*se4[4]), coef1[19]-(1.96*se4[5]),
            coef1[20]-(1.96*se4[6]))

study <- c(rep('Study I', 10), rep('Study II', 10))
party <- rep(c('No', 'No', 'No', 'No', 'No',
               'Yes', 'Yes', 'Yes', 'Yes', 'Yes'), 2)
tr <- rep(c('eAnti-Vax', 'dDemocrat', 'cRepublican',
            'bAnti-Vax Democrat', 'aAnti-Vax Republican'), 4)

df_4 <- as.data.frame(cbind(coef1, hi1, lo1, party, tr, study))
df_4[, 1:3] <- lapply(df_4[, 1:3], as.numeric)
df_4[, 4:6] <- lapply(df_4[, 4:6], as.factor)
df_4$stupid <- str_c(df_4$party, ', ', df_4$study)

# creating figure G.2
plot4 <- ggplot(df_4, aes(shape = party)) +
  theme_bw() +
  geom_vline(xintercept = 0, linetype = 'dashed') +
  geom_pointrange(aes(x = coef1, y = tr, xmin = lo1,
                         xmax = hi1), size = 1,
                     lwd = 1, position = position_dodge(width = 0.5), fill = 'WHITE') +
  scale_shape_manual(name = 'Close Covid Death',
                    labels = c('Yes', 'No'),
                    values = c(15, 1)) +
  scale_y_discrete(name = 'Treatment Condition',
                   labels = c('Anti-Vax Republican', 'Anti-Vax Democrat',
                              'Republican', 'Democrat', 'Anti-Vax')) +
  scale_x_continuous("Estimated Effect on Sympathy for Victim") +
  facet_wrap(~study) +
  labs(title = 'Group Empathy and Sympathy',
       caption = '\nNote: Estimated effects were calculated from two separate models (based on partisan identity). \nError bars represent 95/% confidence intervals. Full results table in Appendix E.') +
  theme(axis.text=element_text(size = 10)) + 
  theme(plot.title = element_text(hjust = 0.5)) + 
  theme(plot.caption = element_text(hjust = 0)) +
  theme(text=element_text(size = 10, family="serif"))

plot4_sl <- ggplot(df_4, aes(shape = party)) +
  theme_bw() +
  geom_vline(xintercept = 0, linetype = 'dashed') +
  geom_pointrange(aes(x = coef1, y = tr, xmin = lo1,
                         xmax = hi1), size = 1,
                     lwd = 1, position = position_dodge(width = 0.5), fill = 'WHITE') +
  scale_shape_manual(name = 'Close Covid Death',
                    labels = c('Yes', 'No'),
                    values = c(15, 1)) +
  scale_y_discrete(name = 'Treatment Condition',
                   labels = c('Anti-Vax Republican', 'Anti-Vax Democrat',
                              'Republican', 'Democrat', 'Anti-Vax')) +
  scale_x_continuous("Estimated Effect on Sympathy for Victim") +
  facet_wrap(~study) +
  labs(title = 'Group Empathy and Sympathy') +
  theme(axis.text=element_text(size = 10)) + 
  theme(plot.title = element_text(hjust = 0.5)) + 
  theme(plot.caption = element_text(hjust = 0)) +
  theme(text=element_text(size = 10, family="serif"))

plot1; plot2; plot3; plot4

# ggsave('plot1.png', plot1, 'png', '~/Dropbox/Wayde/Research/PostDiss_Trauma/CES21/Study1_Covid',
#     width = 190, height = 127, units = 'mm')
# ggsave('plot2.png', plot2, 'png', '~/Dropbox/Wayde/Research/PostDiss_Trauma/CES21/Study1_Covid',
#     width = 190, height = 127, units = 'mm')
# ggsave('plot3.png', plot3, 'png', '~/Dropbox/Wayde/Research/PostDiss_Trauma/CES21/Study1_Covid',
#     width = 190, height = 127, units = 'mm')
# ggsave('plot4.png', plot4, 'png', '~/Dropbox/Wayde/Research/PostDiss_Trauma/CES21/Study1_Covid',
#     width = 190, height = 127, units = 'mm')

# ggsave('plot1_sl.png', plot1_sl, 'png', '~/Dropbox/Wayde/Research/PostDiss_Trauma/CES21/Study1_Covid',
#     width = 190, height = 127, units = 'mm')
# ggsave('plot2_sl.png', plot2_sl, 'png', '~/Dropbox/Wayde/Research/PostDiss_Trauma/CES21/Study1_Covid',
#     width = 190, height = 127, units = 'mm')
# ggsave('plot3_sl.png', plot3_sl, 'png', '~/Dropbox/Wayde/Research/PostDiss_Trauma/CES21/Study1_Covid',
#     width = 190, height = 127, units = 'mm')
# ggsave('plot4_sl.png', plot4_sl, 'png', '~/Dropbox/Wayde/Research/PostDiss_Trauma/CES21/Study1_Covid',
#     width = 190, height = 127, units = 'mm')

```

## Alternative Modeling Approaches

```{r}

df$tr_vax <- ifelse((df$tr1 == 2 | df$tr1 > 4 ), 1, 0)
df$tr_out <- ifelse(df$gop == 1 & (df$tr1 == 3 | df$tr1 == 5), 1,
             ifelse(df$gop == 0 & (df$tr1 == 4 | df$tr1 == 6), 1, 0))
df2$tr_vax <- ifelse((df2$trgrp == 2 | df2$trgrp > 4 ), 1, 0)
df2$tr_out <- ifelse(df2$gop == 1 & (df2$trgrp == 3 | df2$trgrp == 5), 1,
             ifelse(df2$gop == 0 & (df2$trgrp == 4 | df2$trgrp == 6), 1, 0))

m_add_blm <- lm(study1_dv_1 ~ tr_vax*tr_out, data = df) 
summary(m_add_blm)
m_add_sym <- lm(study1_dv2_1 ~ tr_vax*tr_out, data = df) 
summary(m_add_sym)
m_add_blm2 <- lm(study1_dv_1 ~ tr_vax*tr_out, weights = teamweight, data = df2) 
summary(m_add_blm2)
m_add_sym2 <- lm(study1_dv2_1 ~ tr_vax*tr_out, weights = teamweight, data = df2) 
summary(m_add_sym2)

m_add_blm_r <- lm(study1_dv_1 ~ tr_vax*tr_out, data = subset(df, gop == 1)) 
summary(m_add_blm_r)
m_add_sym_r <- lm(study1_dv2_1 ~ tr_vax*tr_out, data = subset(df, gop == 1)) 
summary(m_add_sym_r)
m_add_blm_r2 <- lm(study1_dv_1 ~ tr_vax*tr_out, weights = teamweight, 
                   data = subset(df2, gop == 1)) 
summary(m_add_blm_r2)
m_add_sym_r2 <- lm(study1_dv2_1 ~ tr_vax*tr_out, weights = teamweight,
                   data = subset(df2, gop == 1)) 
summary(m_add_sym_r2)

m_add_blm_d <- lm(study1_dv_1 ~ tr_vax*tr_out, data = subset(df, gop == 0)) 
summary(m_add_blm_d)
m_add_sym_d <- lm(study1_dv2_1 ~ tr_vax*tr_out, data = subset(df, gop == 0)) 
summary(m_add_sym_d)
m_add_blm_d2 <- lm(study1_dv_1 ~ tr_vax*tr_out, weights = teamweight,
                   data = subset(df2, gop == 0)) 
summary(m_add_blm_d2)
m_add_sym_d2 <- lm(study1_dv2_1 ~ tr_vax*tr_out, weights = teamweight,
                   data = subset(df2, gop == 0)) 
summary(m_add_sym_d2)

# creating table D.4
modelsummary(models = list('Blame'  = m_add_blm,
                           'Sympathy' = m_add_sym,
                           'Blame'  = m_add_blm2,
                           'Sympathy' = m_add_sym2),
             output = 'latex',
             stars = c('*' = .05),
             gof_omit = 'AIC|BIC|F|RMSE|Log.Lik.',
             escape = FALSE)

# creating table D.5
modelsummary(models = list('Blame'  = m_add_blm_d,
                           'Sympathy' = m_add_sym_d,
                           'Blame'  = m_add_blm_d2,
                           'Sympathy' = m_add_sym_d2),
             output = 'latex',
             stars = c('*' = .05),
             gof_omit = 'AIC|BIC|F|RMSE|Log.Lik.',
             escape = FALSE)

# creating table D.6
modelsummary(models = list('Blame'  = m_add_blm_r,
                           'Sympathy' = m_add_sym_r,
                           'Blame'  = m_add_blm_r2,
                           'Sympathy' = m_add_sym_r2),
             output = 'latex',
             stars = c('*' = .05),
             gof_omit = 'AIC|BIC|F|RMSE|Log.Lik.',
             escape = FALSE)


```


## Code Graveyard

```{r Covid Study, eval = FALSE}

df1_1 <- group.CI(study1_dv_1 ~ tr1,
                   data = df,
                   ci = 0.90)
df1_1_dem <- group.CI(study1_dv_1 ~ tr1,
                   data = subset(df, gop == 0),
                   ci = 0.90)
df1_1_gop <- group.CI(study1_dv_1 ~ tr1,
                   data = subset(df, gop == 1),
                   ci = 0.90)

df1_2 <- group.CI(study1_dv_2 ~ tr1,
                   data = df,
                   ci = 0.90)
df1_2_dem <- group.CI(study1_dv_2 ~ tr1,
                   data = subset(df, gop == 0),
                   ci = 0.90)
df1_2_gop <- group.CI(study1_dv_2 ~ tr1,
                   data = subset(df, gop == 1),
                   ci = 0.90)

df1_3 <- group.CI(study1_dv_3 ~ tr1,
                   data = df,
                   ci = 0.90)
df1_3_dem <- group.CI(study1_dv_3 ~ tr1,
                   data = subset(df, gop == 0),
                   ci = 0.90)
df1_3_gop <- group.CI(study1_dv_3 ~ tr1,
                   data = subset(df, gop == 1),
                   ci = 0.90)

df1_sym <- group.CI(study1_dv2_1 ~ tr1,
                   data = df,
                   ci = 0.90)
df1_sym_dem <- group.CI(study1_dv2_1 ~ tr1,
                   data = subset(df, gop == 0),
                   ci = 0.90)
df1_sym_gop <- group.CI(study1_dv2_1 ~ tr1,
                   data = subset(df, gop == 1),
                   ci = 0.90)

names(df1_1) <- c('group', 'hi', 'mean', 'lo')
names(df1_1_dem) <- c('group', 'hi', 'mean', 'lo')
names(df1_1_gop) <- c('group', 'hi', 'mean', 'lo')

names(df1_2) <- c('group', 'hi', 'mean', 'lo')
names(df1_2_dem) <- c('group', 'hi', 'mean', 'lo')
names(df1_2_gop) <- c('group', 'hi', 'mean', 'lo')

names(df1_3) <- c('group', 'hi', 'mean', 'lo')
names(df1_3_dem) <- c('group', 'hi', 'mean', 'lo')
names(df1_3_gop) <- c('group', 'hi', 'mean', 'lo')

names(df1_sym) <- c('group', 'hi', 'mean', 'lo')
names(df1_sym_dem) <- c('group', 'hi', 'mean', 'lo')
names(df1_sym_gop) <- c('group', 'hi', 'mean', 'lo')

df1_1$dv <- rep('Victim', 6)
df1_1_dem$dv <- rep('Victim', 6)
df1_1_gop$dv <- rep('Victim', 6)

df1_2$dv <- rep('Govt', 6)
df1_2_dem$dv <- rep('Govt', 6)
df1_2_gop$dv <- rep('Govt', 6)

df1_3$dv <- rep('Activists', 6)
df1_3_dem$dv <- rep('Activists', 6)
df1_3_gop$dv <- rep('Activists', 6)

df1_sym$dv <- rep('Sympathy', 6)
df1_sym_dem$dv <- rep('Sympathy', 6)
df1_sym_gop$dv <- rep('Sympathy', 6)

dat1 <- bind_rows(df1_1,df1_2, df1_3, df1_sym)
dat1_dem <- bind_rows(df1_1_dem, df1_2_dem, df1_3_dem, df1_sym_dem)
dat1_gop <- bind_rows(df1_1_gop, df1_2_gop, df1_3_gop, df1_sym_gop)

```

```{r echo=TRUE, eval = FALSE}

## Plotting

cbPalette <- c("#CC79A7", "#009E73", "#F0E442", "#0072B2", "#D55E00")

plot1.1 <- ggplot(data = dat1, aes(x = group, y = mean, fill = dv)) +
  theme_bw() + scale_fill_manual(values = cbPalette) +
  geom_bar(stat="identity", position=position_dodge()) +
  geom_errorbar(aes(ymin = lo, ymax = hi),
                size = 0.5, width = 0.1,
                position = position_dodge(0.9)) + 
  scale_x_discrete(name = '\nTreatment Condition', 
                   labels = c('Control', 'Anti-Vax', 'Democrat', 'Republican', 'Anti-Vax Dem', 'Anti-Vax Rep')) +
  coord_cartesian(ylim=c(20, 85)) +
  ylab('Average Blame/Sympathy') +
  labs(title = "Effect of Party and Anti-Vax on Blame Attribution") +
  theme(plot.title = element_text(hjust = 0.5)) +
  theme(legend.position="bottom") +
  theme(text=element_text(family="serif")) +
  guides(fill=guide_legend(title="Blame/Sympathy"))
plot1.1

plot1.dem <- ggplot(data = dat1_dem, aes(x = group, y = mean, fill = dv)) +
  theme_bw() + scale_fill_manual(values = cbPalette) +
  geom_bar(stat="identity", position=position_dodge()) +
  geom_errorbar(aes(ymin = lo, ymax = hi),
                size = 0.5, width = 0.1,
                position = position_dodge(0.9)) + 
  scale_x_discrete(name = '\nTreatment Condition', 
                   labels = c('Control', 'Anti-Vax', 'Democrat', 'Republican', 'Anti-Vax Dem', 'Anti-Vax Rep')) +
  coord_cartesian(ylim=c(20, 85)) +
  ylab('Average Blame/Sympathy') +
  labs(title = "Effect of Party and Anti-Vax on Blame Attribution, Democrats only") +
  theme(plot.title = element_text(hjust = 0.5)) +
  theme(legend.position="bottom") +
  theme(text=element_text(family="serif")) +
  guides(fill=guide_legend(title="Blame/Sympathy"))
plot1.dem

plot1.gop <- ggplot(data = dat1_gop, aes(x = group, y = mean, fill = dv)) +
  theme_bw() + scale_fill_manual(values = cbPalette) +
  geom_bar(stat="identity", position=position_dodge()) +
  geom_errorbar(aes(ymin = lo, ymax = hi),
                size = 0.5, width = 0.1,
                position = position_dodge(0.9)) + 
  scale_x_discrete(name = '\nTreatment Condition', 
                   labels = c('Control', 'Anti-Vax', 'Democrat', 'Republican', 'Anti-Vax Dem', 'Anti-Vax Rep')) +
  coord_cartesian(ylim=c(20, 85)) +
  ylab('Average Blame/Sympathy') +
  labs(title = "Effect of Party and Anti-Vax on Blame Attribution, Republicans only") +
  theme(plot.title = element_text(hjust = 0.5)) +
  theme(legend.position="bottom") +
  theme(text=element_text(family="serif")) +
  guides(fill=guide_legend(title="Blame/Sympathy"))
plot1.gop

```

